Estimating Access to Rabies Post-Exposure Prophylaxis in Madagascar
INTRODUCTION
Deaths due to canine mediated rabies, estimated to cause approximately 60,000 human deaths anually, can be prevented through prompt administration of post-exposure prophylaxis. However, access to this intervention is highly limited in areas where the disease is endemic, often due to inaccessibility to health centres which provide PEP (cite GAVI paper or Nandini’s paper of geographic availability of PEP). Data on true rabies exposures in humans and incidence in animals is also lacking in most of these countries, with the most commonly available data being numbers of bite victims reporting to health facilities.
Here–discuss geographic variation in access to vaccination and care–how this shapes mortality for other diseases.
The majority of rabies burden studies use these data to estimate burden of rabies from probability decision tree frameworks, often with the key assumption that overall reported bite incidence (i.e. both bites due to non-rabid and rabid animals) are proportional to rabies incidence (i.e. the more bites reported in a location, the higher the incidence of rabies exposures there) and that reporting is uniform across space. While at the national level these estimates may be accurate, at the sub-national level, this framework will likely underestimate rabies deaths in places with low reporting and overestimates rabies deaths in places with high reporting of bites. The most recent estimation of burden and the impact of PEP used another approach, using transmission dynamic models as the backbone to predict incidence based on the level of vaccination coverage and size of the dog population at the national level. Using transmission dynamic models to estimate incidence could improve upon previous studies which may underestimate rabies burden in areas with low reporting.
In Madagascar, Institut Pasteur de Madagascar (IPM) provides PEP at no-cost to patients at 31 clinics across the country. PEP is not available at any other public clinics or through the private sector. In addition, there is limited control in dog populations and the disease is endemic throughout the country. Due to the spatially restricted nature of PEP and the endemic circulation of canine rabies, geographic access is likely to be a major driver of disease burden within the country. To get spatial estimates of disease burden in this context, we flip the standard decision tree and make the assumption that reported bite incidence reflects access to PEP rather than differences in rabies incidence, using travel times to clinics as a predictor of reported bites. Then using a range of rabies incidence given endemic transmission with no mass dog vaccination (GAVI paper), we generate sub-national estimates of rabies burden in an adapted decision tree framework. Finally, using this same model pipeline, we explore the impacts of geographically expanding access to PEP in Madagascar on reducing human rabies deaths.
METHODS
GIS Data
We used the global friction surface for 2015 generated by the Malaria Atlas Project ( https://map.ox.ac.uk/research-project/accessibility_to_cities/, Weiss et al. 2015,) and GPS points of clinics to estimate the travel time to the nearest ARMC for the country of Madagascar at a 1 x 1 km scale. We then calculated a weighted average of travel times by human population to the commune level, using administrative shapefiles available trhough the UN Office for the Coordination of Humanitarian Affairs. Human population estimates were taken from the 2015 UN adjusted population projections from World Pop (www.worldpop.org, Linaird et al. 2012) and also aggreagated to the commune level.
Bite patient data
We used a database of bite patient forms submitted to IPM from ARMC across the country between 2014 - 2017. These were individual patient data forms that were submitted as frequently as monthly to annually by clinics. Of the 31 existing ARMC, 10 submitted fewer than 10 forms over the four years. Two clinics, the IPM ARMC and the Fianarantsoa ARMC had separate databases which were not available at the time of analysis. Overall, we had data from 19 clinics across the country (Fig 1A). These data include the administrative district of the bite patient’s address and the date of reporting. We also had 27 months of data from the Moramanga District that were resolved to the commune level (the administrative level below district).
For most districts, the majority of bites were reported to the closest clinic as estimated by our travel time metric (Fig 1A, Fig S1). Therefore, we defined the catchment area for each clinic as all districts for which the clinic was the closest ARMC. We excluded any districts from catchments for which the clinic did not submit any forms (n = 12, districts in grey in Fig 1A).
Figure 1A
As even for clinics which submitted data, there was substantial undersubmission of forms, we estimated clinic level reporting as the proportion of days on which forms were submitted, excluding any periods for which there were no forms submitted for 15 consecutive days (Figure S2 shows the time series of form submission for each clinic and periods of time excluded by the 15 day threshold). Our estimate of reporting did vary based on our assumption of the threshold number of consecutive days (Figure S3), so we did look at the sensitivity of model parameter estimates to changing this threshold. To estimate the average annual bites reported for each district, we further excluded data from any years for which there was less than 25% reporting at the clinic level.
Figure 2A
Figure 1B
Figure 3B
Figure 3A
Previous work in the Moramanga District showed that low risk contacts with probable or confirmed rabies cases make up approximately 20% of patients reporting to ARMC (Rajeev et al. 2018). Generally, contacts present as clustered cases, so we excluded patients reporting on any dates that had greater than 3 standard deviations above the mean number of patients reporting per day (Figure 3A). We validated this method using the Moramanga ARMC data for which we had details on the type of exposure, and found that setting the threshold to 3 standard deviations (SDs) resulted in approximately 50% of known contacts excluded, with only 2% of non-contacts excluded (Figure 3B). For the national data for which a subset of patient forms were explicitly noted to be contacts, we found that our exclusion criteria of 3 SDs identified 68.28 of known contacts. We further excluded these known contacts which were not identified based on the daily distribution of patients, resulting in the exclusion of approximately 7.18% of patient records from the national data. We also compared this method to assuming that 20% of total bites at the district level were contacts.
Figure 4
After excluding contacts and correcting for undersubmission of forms, our final dataset consisted of estimates of average bite incidence for 71 districts from 19 catchments (Figure 4).
Model of reported bites as a function of travel time
While the national bite patient data is available at the district level, travel times and therefore reporting, likely vary significantly at the sub-district level. In order to translate the impacts of differences in access at sub-district scales to the magnitude of reported bites at the district scale, we modeled bites at the district level as the sum of incidence at the commune level. Incidence at the commune level is then a function of travel times to the closest ARMC. Specifically, we modeled bites as follows:
\[ \mu_{d} = \sum \limits_{j=1}^jexp(\beta_{t}T_j + \beta_0)\times pop_j \] where \(\mu_{d}\) is the mean number of bites in district, which is the sum of bites at the commune level given commune level travel times, \(T_j\). We then estimate the likelihood of observing the bites at the district level where bites are a poisson distribution around the mean \(\mu_{d}\), given \(\beta_t\), the effect of travel times of reported bites and \(\beta_0\), the model intercept.
As we had data available on reported bites at the commune level from the Moramanga ARMC (from Rajeev et al. 2018), we modeled observed commune bites in the same framework, but un-aggregated, where:
\[ \mu_{j} = exp(\beta_tT_j + \beta_0)\times pop_j \]
where \(\mu_{j}\) = the mean number of bites in commune \(j\) and the observed bites at the commune level follow a poisson distribution around the mean \(\mu_d\). We only included communes which were designated to be within the catchment for the clinic. We compared our estimates of \(\beta_t\) (i.e. the impact of travel time on incidence) and \(\beta_0\) (the intercept) for our district data at the national level and the commune level data from the Moramanga ARMC.
Finally, we compared these models to a model of bites at the district level as the function of travel times averaged to the district level: \[ \mu_{d} = exp(\beta_{t}T_d + \beta_0)\times pop_d \] where \(pop_d\), \(T_d\), and \(mu_d\) are respectively human population size, weighted travel times, and mean number of bites in district \(d\). As travel times are correlated with human population size (Figure Z), we also compared how well bites were predicted by human population size alone for these different models, as a test of whether the observed patterns could be predicted by bite incidence scaling with population size. We also looked at how well distance from the closest ARMC (rather than travel times) predicted bites as an alternative proxy for access.
Estimation of burden and reporting
We used our model to predict average annual bite incidence for all 114 districts in Madagascar, and estimated average reporting of rabid exposures and deaths due to rabies given this and assumptions about rabies exposures.
We calculated the expected reporting of rabid exposures (\(\overline\rho\)) given bite incidence as predicted by our model(\(\mu\)) as: \[
\overline\rho = \frac{\mu \times p_{rabid}}{R}
\] or the fraction of incidence due to rabid exposures (\(\mu \times p_{rabid}\)) divided by the total rabies exposure incidence (\(r\)) for a range of estimated rabies incidence and \(p_{rabid}\). We look at the range of \(p_{rabid}\) reported in Rajeev et al. 2018 for data from the Moramanga District (0.2 - 0.6). where the proportion of reported bites that are rabies exposures (\(p_{rabid}\)) are defined as: \[
p_{rabid}=
\begin{cases}
x, & \text{if}\ \frac{R \times \rho_{max}}{B} > x \\
\frac{R \times \rho_{max}}{B}, & \text{otherwise}
\end{cases}
\]
such that rabid reported bites (i.e. \(p_{rabid} \times B_i\)) cannot exceed the expected number of human exposures given maximum reporting (i.e. \(R_i \times \rho_{max}\)). \(\rho_{max}\) is taken from the Moramanga ARMC data for Moramanga Ville, the commune with the ARMC (i.e. the area with the minimum travel time in the district, on average xx minutes).
To generate \(R\), the rabies incidence in dogs in the absence of any vaccination, \(r\), is multiplied by the estimated dog population in the commune (\(D_i\)) and the exposure rate per rabid dogs (\(p_{exp}\) = 0.39 persons exposed per rabid dog)(Hampson et al. 2018). We use a human:dog ratio (HDR) of 5 to generate our maximum expected incidence and an HDR of 25 for our minimum expected incidence. As there is little data on dog populations in Madagascar, this range of HDRs encompasses a wide range observed across Africa (cite!).
To estimate the average number of deaths for each commune, we extended the above framework into a stochastic framework as follows: \[
deaths_i = (R_i - p_{rabid}B_i) \times p_{death}
\] where \(R_i\) is drawn from a uniform distribution between the minimum and maximum expected number of human exposures for each location and \(B_i\), the number of reported bites, is drawn from a poisson distribution with the mean predicted number of bites from the travel time model. We constrain \(p_rabid\) as in Equation N, and we assume that all rabies exposures reported to an ARMC receive and complete PEP, and PEP is completely effective at preventing death due to rabies. The probability of death in the absence of PEP is taken from cite GAVI/Joel paper.
Estimating the impact of expanding PEP Access (TO DO)
We use this framework to compare three scenarios of PEP provisioning in Madagascar:
- The baseline with the current clinic locations (n = 31)
- Expansion of ARMC to one clinic per district (n = 114)
- Expansion of ARMC to two clinics per district, choosing the clinic which minimizes the proportion of people living greater than 2 hours away from a clinic.
We use data on the location of clinics provided by IPM to regenerate travel times to the nearest ARMC given expansion as per scenario 2 and 3. We then predict the expected bite incidence from the model given these new travel times and compare the relative decreases in burden for the three scenarios.
Strategically expanding PEP access (TO DO)
Given limited resources and capacity of clinics to provision PEP, we developed a framework to look at the incremental benefit of additional ARMC. Starting with the current locations (Scenario 1), we added one clinic and recalculated the proportion of people living < 3 hours away from a clinic for the country. We added the clinic which minimized this metric, and then repeated the process iteratively, ranking clinics and adding the top clinic in each interation sequentially. We calculate burden for the given clinic locations and look at the incremental reduction in burden as each clinic is added.
Sensitivity analyses for burden estimates and scenario analyses (TO DO)
RESULTS
Models of bites as a function of travel times
Figure 5A
| model | par | values | upper_CI | lower_CI |
|---|---|---|---|---|
| Mada communes | B_ttimes | -0.0098051 | -0.0096175 | -0.0099809 |
| Moramanga communes | B_ttimes | -0.0125278 | -0.0113804 | -0.0137029 |
| Mada districts | B_ttimes | -0.0082565 | -0.0081410 | -0.0083732 |
| Mada communes | intercept | -4.9394093 | -4.9231004 | -4.9569247 |
| Moramanga communes | intercept | -5.2735226 | -5.1566258 | -5.3818631 |
| Mada districts | intercept | -5.3488643 | -5.3479843 | -5.3497870 |
Figure 5B
We estimated similar parameter values from our commune-level data from the Moramanga ARMC and the district level data from 19 clinics across the country (Table 1, Figure 5A), with reported bite incidence decreasing with travel times to the ARMC. All of the models produced reasonable fits to the data (Figure 5B), however, there was some variation in bite incidence that was not captured by the model.
Model validation: TO DO
Estimation of burden and reporting
Figure 6
Generally, estimated reporting of rabies exposures decayed with travel times given model predicted bite incidence and a range of rabies incidence and \(p_{rabid}\) (Figure 6). Given our model assumptions, reporting was estimated at the maximum of 0.98 for travel times under 1 hour given the maximum estimated rabies exposure incidence and the minimum estimate of \(p_{rabid}\) (the lower range of reporting probabilities), and travel times under 5.5 hours given the minimum estimated rabies exposure incidence and the maximum estimate of \(p_{rabid}\) (the upper range of reporting probabilities).
Figure 7
Figure 7
Figure 7
Figure 7
When we estimate burden of deaths stochastically within this range of incidence and given our high and low estimates of proportion of reported bites that are rabid, we see that burden of deaths also decreases with travel times (Figure 7, results presented aggregated at the district level). Overall, we estimate average annual deaths between 383 - 705. Also estimate deaths averted here!.
When we compare our burden estimates at the district vs. the commune level (summed to district), we see that while overall, the estimates are very similar (Fig NA), at low travel times, calculating burden at the district level results in an assumption of maximum reporting for the whole district, which assumes very low burden
Sensitivity analyses
DISCUSSION
Key findings
Strengths and Limitations
- Not accounting for clinic functioning
- Not explicit data on reporting or rabies incidence!
- Other factors that drive reporting
- While we use a signficiant amount of assumptions, these are the same set of assumptions that underlie current rabies estimates, we’re just improving them a bit… and ours are a likely a little bit more in line with reality @ a sub-national level?